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Abstract 

The uncoupled Continuous Time Random Walk {CTRW) in one 
space-dimension and under power law regime is splitted into three 
distinct random walks: (rw\), a random walk along the line of 
natural time, happening in operational time; (rw-z), a random walk 
along the line of space, happening in operational time; (rws), the 
inversion of (rw\), namely a random walk along the line of operational 
time, happening in natural time. Via the general integral equation 
of CTRW and appropriate rescaling, the transition to the diffusion 
limit is carried out for each of these three random walks. Combining 
the limits of (rw\) and (ru^) we get the method of parametric 
subordination for generating particle paths, whereas combination of 
(rw2) and (rw^) yields the subordination integral for the sojourn 
probability density in space - time fractional diffusion. 
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1. Introduction 



In recent decades the notion of continuous time random walk has become 
popular as a conceptual and computational model for diffusion processes. 
Its relation to continuous diffusion processes is an asymptotic one, either 
by considering spatial and temporal refinement or by looking for large-time 
wide-space behaviour (both views being equivalent to change the units of 
measurement in time and space). Often such walks are considered for power 
law behaviour of waiting times and jumps, and from the behaviour in the 
Laplace-Fourier domain near zero the relevant limiting diffusion equations 
are deduced heuristically. There are several methods of making precise this 
transition to the limit, the fundamental role being played by the domains of 
attraction of stable probability laws (in time and space). It is desirable 
to have presentations of this transition without much theory of abstract 
probability and stochastic processes, and this can be done in a transparent 
manner for the derivation of the limiting fractional diffusion equations. The 
essential idea of our approach presented here is to use the integral equation of 
continuous time random walk separately for three random walks obtained by 
splitting the original walk, and then synthesizing the results. We so obtain 
our method of "parametric subordination" as well as what usually is called 
"subordination" in fractional diffusion. 

2. The Space-Time Fractional Diffusion 

We begin by considering the Cauchy problem for the (spatially one- 
dimensional) space-time fractional diffusion (STFD) equation. 

t D?u(x,t) = x Dffu{x,t), u(x,0) = 6(x), x e M, t>0, (2.1) 

where {a , 6 , (3} are real parameters restricted to the ranges 

0<«<2, |0| < min{a,2 - a} , < (3 < 1 . (2.2) 

Here t D@ denotes the Caputo fractional derivative of order j3, acting on the 
time variable t, and x Dq denotes the Riesz-Feller fractional derivative of 
order a and skewness 6, acting on the space variable x. Let us note that the 
solution u(x,t) of the Cauchy problem (2.1), known as the Green function 
or fundamental solution of the space-time fractional diffusion equation, is a 
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probability density in the spatial variable x, evolving in time t. In the case 
a = 2 and /3 — 1 we recover the standard diffusion equation for which the 
fundamental solution is the Gaussian density with variance a 2 = It. 

Writing, with Re[s] > cr , k G 1R, the transforms of Laplace and Fourier as 

POO 

C{f(t);s} = f(s):= / e- st f(t)dt, 

Jo 

/+oo . 
e lKX g(x) dx, 
-oo 

we have the corresponding transforms of t D& fit) and x D^g(x) as 

C{ t DPf(t)} = sPf(s)- S P- 1 f(0), (2.3) 

F{ x D?g(x)} = -\K\ a i e si S n K g{n). (2.4) 

Notice that ^ s ig nK = exp [i (signre) 9ir/2\. For the mathematical details 
the interested reader is referred to Gorenflo and Mainardi (1997), Podlubny 
(1999), Butzer and Westphal (2000), Kilbas, Srivastava and Trujillo 
(2006), Hilfer (2008), Tomovski, Hilfer and Srivastava (2010) on fractional 
derivatives, and to Samko, Kilbas and Marichev (1993), Rubin (1996), on 
the Feller potentials. 

For our purposes let us here confine ourselves to recall the representation 
in the Laplace- Fourier domain of the (fundamental) solution of (2.1) as it 
results from the application of the transforms of Laplace and Fourier. Using 
5(k) e 1 we have from (2.1) 

s $ s) - s f 3 - 1 = - 1 k | a % 6 si S n K s) , 

hence 

u(k, s) = —5 q—. . (2.5) 

sP + \k\ a i" sl S n K 

For explicit expressions and plots of the fundamental solution of (2.1) in the 
space-time domain we refer the reader to Mainardi, Luchko and Pagnini 
(2001). There, starting from the fact that the Fourier transform u{K,t) 
can be written as a Mittag-Leffler function with complex argument, the 
authors have derived a Mellin-Barnes integral representation of u(x,t) with 
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which they have proved the non- negativity of the solution for values of the 
parameters {a, 9, (3} in the range (2.2) and analyzed the evolution in time 
of its moments. In particular for {0 < a < 2, (3 = 1} we obtain the stable 
densities of order a and skewness 9. The representation of u(x,t) in terms 
of Fox if-functions can be found in Mainardi, Pagnini and Saxena (2005), 
see also Chapter 6 in the book by Mathai, Saxena and Haubold (2010). We 
note, however, that the solution of the STFD Equation (2.1) and its variants 
has been investigated by several authors; let us only mention some of them, 
Barkai (2001), Meerschaert et al. (2002 and 2004), Metzler and Klafter 
(2004), where the connection with the CTRW was also pointed out. 



3. The Continuous-Time Random Walk 



Starting in the Sixties and Seventies of the past century the concept of 
continuous time random walk, CTRW, became popular in physics as a rather 
general (microscopic) model for diffusion processes, see the monograph by 
Weiss (1994). Mathematically, a CTRW is a renewal process with rewards 
or a random-walk subordinated to a renewal process, and has been treated as 
such by Cox (1967). It is generated by a sequence of independent identically 
distributed (iid) positive waiting times T\,T 2 ,T 3 , . . . , each having the same 
probability distribution $(£),£> 0, of waiting times T , and a sequence of iid 
random jumps X ± , X 2 , X 3 , . . . , each having the same probability distribution 
function W(x), x e IR, of jumps X. These distributions are assumed to be 
independent of each other. Allowing generalized functions (interpretable as 
measures), we have corresponding probability densities <f)(t) = $'(£) and 
w(x) = W'(x) that we will use for ease of notation. Setting t = 0, 

t n = Ti + T 2 + h T n for n G IN, and x = 0, x n = X 1 + X 2 + h X n , 

x{t) = x n for t n < t < t n+ i, we get a (microscopic) model of a diffusion 
process. A wandering particle starts in xq = and makes a jump X n at 
each instant t n . Natural probabilistic reasoning then leads us to the integral 
equation of continuous time random walk for the probability density p(x, t) 
of the particle being in position x at instant t > 0: 



p(x,t)=5(x)^(t) + jy t -t') j 



w(x — x') p(x', t') dx' 



dt! , (3.1) 
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Here the survival probability 

/oo 
0(0 dt' (3.2) 

denotes the probability that at instant t the particle still is sitting in its 
initial position xq = 0. Eq. (3.1) is also known as (semi -) Markov renewal 
equation. Germano et al. (2009), considering the more general situation of 
coupled CTRW, have given a detailed stochastic proof of an integral equation 
which by modification and specialization reduces to (3.1). 

By using the Fourier and Laplace transforms we arrive, via S(k) = 1 and the 
convolution theorems, in the transform domain at the equation 

p(k, s) = \&(s) + w(k) (p(s) p(k, s) , (3.3) 

which, by \f(s) = (1 — 0(s))/s implies the Montroll- Weiss equation, see e.g 
Weiss (1994), 

p(«,a) = — ~~— ■ (3.4) 

s 1 — w(k) <p{s) 

Because of \w(k)\ < 1, \4>(s)\ < 1 for k ^ 0, s ^ 0, we can expand 

oo 

= [?(s) w(k)Y , (3.5) 

n=0 

and promptly obtain, see e.g. Cox (1967) and Weiss (1994), 

oo 

p(x, t) = v n (t) w n {x) — oo < x < +oo, t > . (3.6) 

ra=0 

Here the functions v n = ^ * 0* n and u> n = u>* n are obtained by iterated 
convolutions in time t and in space x, respectively; in particular we have: 

v (t) = (* * S)(t) = «i(t) = (* * 0)(O, ™o(z) = <5(x), fiW = ™(z). 

(3-7) 

The representation (3.6) can be found without the detour over (3.4) by 
direct probabilistic treatment. It exhibits the CTRW as a subordination 
of a random walk to a renewal process. 
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Note that in the special case <fi(t) = mexp (—mt), m > 0, the equation (3.1) 
describes the compound Poisson process. It reduces after some manipulations 
to the Kolmogorov- Feller equation 



d_ 

dt 1 



-p(x,t) = —mp(x,t) + m / w(x — x')p(x', t) dx' . (3.8) 
From (3.6) we then obtain the well-known series representation 

p( x ,t)=e~ mt f2 { -^w n (x). (3.9) 



n=0 



4. Power Law Assumptions and Well- Scaled 
Passage to the Diffusion Limit 

We will pass to the diffusion limit of CTRW under power law assumptions 
for waiting times and jumps given by asymptotic relations 

l-0(s)~As /3 , A>0, s^0 + , (4.1) 

1 - w(k) ~ i2\K\ a i ds[ S nK , fl>0, k^O, (4.2) 

where /3, a and 9 are restricted as in (2.2). See, e.g. Gorenflo and Mainardi 
(2003), Scalas, Gorenflo and Mainardi (2004), Gorenflo (2010) for their 
physical meaning, obtainable via Tauber-Karamata theorems on aymptotics. 
We rescale space and time by multiplying the jumps Xk and waiting times 
Tk by positive factors h and r (intended to be sent to zero), so obtaining a 
rescaled walk 

n n 

x n (h) = ^2hX k , t„(r) = ^rT fc . 

k=l k=l 

Physically, this amounts to changing the units of measurement from 1 to 1/h 
in space x, from 1 to 1/r in time t. 

For the reduced jumps and waiting times we get 

w h (x) = w(x/h)/h , <t> T {t) = (j)(t/r)/T , (4.3) 
w h (K)=w(hK),, T (s) = 0(rs), (4.4) 
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and the Montroll- Weiss equation (3.3) transforms into 



g h>T («, 8 ) = j-^ - 1 „ = _J__ , (4 . 5) 

s 1 — Wh{n) 4> T (s) s 1 — w(hK) (f>(rs) 

Using our power law assumptions (4.1), (4.2) we find, for h — > 0, r — >■ 0, the 
asymptotics 



-i 



p hT (K, s) ~ g-^ = g-^ (4.6) 

//(/i|«;|)°^ sl g n « + A(rs)^ p(/i,r) i d sl S n K + 

with 

p{h,T) 



Adopting the scaling relation 



P(h,r) = ^ = l, (4.7) 



we arrive for h — > at the limit 



p (/c, s) = q—. =: u(k, s) 

which we identify as the Fourier-Laplace solution to the Cauchy problem (2.1) 
for the STFD equation. By inverting the Laplace transform and then for the 
Fourier transform invoking the continuity theorem of probability theory we 
get weak convergence to the fundamental solution of (2.1). 

Remark. Because Mittag-Leffler waiting times distributed according to 
ty(t) = Ep(-tP) play a distinguished role among power law waiting times and 
lead to time-fractional CTRW, in fact to the time-fractional Kolmogorov- 
Feller equation, see e.g. Hilfer (1995), Hilfer and Anton (1995), Gorenflo and 
Mainardi (2008) and Gorenflo (2010), it is natural to use them (properly 
scaled) for approximately simulating particle paths of space-time fractional 
diffusion processes. Fulger, Scalas and Germano (2008) have done this, 
thereby using exact methods for generating the required temporal Mittag- 
Leffler deviates and spatial stable deviates, giving references for these 
methods. 
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5. Splitting into Three Random Walks 



We analyze the anatomy of the CTRW of Section 3 by splitting it into 
three random walks (rwi), (ra 2 ), (rws) that we will treat in analogy to 
our method of well-scaled transition to the diffusion limit. We will work in 
three coordinates: the semiaxis t > of natural time, the semi-axis t* > 
of operational time, the axis — oo < x < +oo of space. In addition to 
the Laplace transform £{/(£); s} = f(s) with respect to t we will need the 
Laplace transform £{f(t*)] s*} = /(s*) with respect to i*. 

For convenience let us denote by (rw) = (CTRW) the random walk of 
Section 3, 

t = 0, t n = T 1 +T 2 + --- + T n ; x Q = 0, x n = X 1 + X 2 + ■ ■ ■ + X n . (5.1) 

This is a random walk along the space line x G 1R happening in physical 
time t > 0. Here (pit) is the waiting time density, w(x) is the jump density, 
and p(x, t) is the sojourn density in point x at instant t. 
By (rwi) we denote the basic renewal process 

Un = n (n = o ? i, 2,...); t = 0, t n = T 1 + T 2 + --- + T n . (5.2) 

We consider it as a random walk along natural time t > 0, happening in 
operational time > 0. Here t* plays the role of a "pseudo time", £ plays 
the role of "pseudo space". The waiting time (in t*) is constant (fixed) =1, 
so we have the "waiting time" density = 8(t* — 1). The jump instants 
are the integers t* n — n > 0, and we set i*o = 0. The "jump density" (in 
t) is 4>(t). We denote the "sojourn density" in "point" t at "instant" by 
r(t,t*). Note that in this sojourn density r the first argument is the pseudo 
space t, whereas the second argument is pseudo time t*. 

By (rw 2 ) we denote the random walk 

t*, n = n(n = 0, 1,2, ...) ; x = , x n = X 1 + X 2 H \- X n , (5.3) 

which we consider as a random walk along natural space x G 1R , happening 
in operational time t* > which (again) plays the role of "pseudo time". 
But x here is the usual spatial variable. The waiting time density in t* is 
(again) v(t*) = 6(t* — l), the jump density is w(x), and we denote the sojourn 
density in point x at "instant" t* by £(x,t*). 
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By (ru> 3 ) we denote the random walk 

f = 0, t n = T 1 +T 2 + --- + T n ; t m = n{n = 0,1, 2,...), (5.4) 

along the semi-axis t* > happening in physical time t > 0. The walk is 
just the inversion of (rwi). The variable t now plays the role of time, the 
waiting time density is <f>(t). But t* plays the role of "pseudo space", the 
jump density being v(t*) = — 1). The sojourn density in "point" t* at 
instant t is denoted by q(t*,t) where the pseudo space t* is the first argument 
and time t is the second argument. 



6. Asymptotic Analysis of the Three Walks 

In (rwi) and (ra 3 ) we have walks in the positive pseudo-space direction. For 
such walks it is convenient to modify the Montroll- Weiss formula by using 
the Laplace transform in place of the Fourier transform. If in the general 
(CTRW) the jump density w(x) > for all x G 1R then with 



g(K) := g(x)e KX dx, re > re , (6.1) 
Jo 

we can replace the original Montroll- Weiss equation (3.4) by 

~/ \ 1 — <t>(s) 1 

p(«,a) = , . ~. . • (6-2) 

s 1 — W{K) <f>(s) 

Here for (rwi) and (ru> 3 ) the specifically relevant replacements of x and t, 
re and s must be made, analogously for (ru^) in (3.3) and (3.4). Analogous 
to Section 4 we will carry out well-scaled transitions to the diffusion limit, 
observing (4.1) and (4.2) and an additional asymptotic relation for the 
Laplace transform of v(t#) = 5(t* — 1). We have 

v(s*) = e~ s * — 1 - s* + 0(3*) , for s* — >• . 

6.1 Analysis of (rwi). 

In place of (6.2) we have 

—/ \ 1 — v(s*) 1 

r(s,s,) = ^ = . (6.3) 

s * 1 — (f>(s) v(s m ) 
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We rescale in t* direction by a factor a > 0, so replacing the constant waiting 
time 1 by a. We replace the "jumps" T by rT with r > 0. So we obtain, for 
t ->■ 0, (J ->■ 0: 

x , , 1 - ^(crs*) 1 



s * 1 — 0(rs) u(<7S*) 

O" 1 



Ar^ + as, M s /3 + s 

a 

By introducing the scaling relation 

Ar^ 



<7 

we arrive for r — >■ at the limit 

r ,o(s, s*) 

whence 



= 1 , (6-4) 



1 



s/ 3 + s* ' 



corresponding to 

r ,o(t,Q = t; 1 ^L/(t; 1 ^t) , (6.5) 

which is a probability density in t evolving with t*. Here LZ (t) denotes the 
right-extremal stable density whose Laplace transform is exp (— s 13 ). 

6.2 Analysis of (ru^). 

In analogy to (3.4) we find 

If \ 1 ~~ ^( s *) 1 fa a\ 

S* 1 —w{k)v(s*) 

and by rescaling (replacing "waiting times" 1 by a > 0, jumps X by hX) 
1 - v(c7S*) 1 



s* 1 — w{hK) v(as*) 

a 1 



ji/i«|K|«i0 si S n K + as, \ K \ a i e si g n K + as, 
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The scaling relation needed here is 

^=1. (6.7) 
o 

and with it we get the limit 



s , + sign*' 
f ,o(«,**)=exp (-q«f ^S™) , 



hence 

£o,o(z, i„) = t; 1 /" (t; 1 /" rr) = / a , e (x, t.) , (6.8) 

which is a probability density in x evolving with t*. Here L a (x) is the 
stable density with Fourier transform exp ^— \n\ a $ s *& n K j and f a> e(x,t) is 
the solution u(x,t) of the Cauchy problem (2.1) for (3 = 1. 

6.3 Analysis of (ru^). 

In place of (6.2) we get 

*A ■ (6.9) 

s l-u(s,)0(s) 

Replacement of waiting times T by rT, "jumps" 1 by <r yields 

1 - 0(rs) 1 



9«t,t(s*,s) 



s 1 - u(<7s*) 0(rs) / 6 1Q \ 



Using the same scaling relation as for (rwi), namely Xr^/a = 1, we get the 
limit 



s* + s 

implying 



go,o( S *' S ) = ~ I „/} » ( 6 - U ) 



5b 1 o(**>s) = ^" 1 e _t * a , (6-12) 
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so that 

5o,o(**>*) = tJ 1 '* 1 r 0fl (t,Q , (6.13) 

with r 0) o as in (6.5) and the operator t J 7 , 7 > of Riemann-Liouville 
fractional integration from to t whose Laplace symbol is s -7 . Note that 
Qo,o{t*,t) is a probability density in t* evolving with t. 

Remark. We also recognize that qo,o(t*,t) can be expressed in terms of the 
so-called M- Wright function, see e.g. the recent book by Mainardi (2010), 
Appendix F, Eqs (F.13), (F.24) and (F.51)-(F.53). In fact, inverting (6.12), 
we get, in view of (F.51) and (F.52), 

q Q ,o{U,t)=r p Mp{U/t p ). (6.14) 

To the same result we arrive inverting (6.11) from s to t, getting 

q 0fi (s*,t)=Ep(-s*t?), (6.15) 

where Ep is the Mittag-Lefrler function. Then, inverting this Laplace 
transform from s* to t* we get, in view of (F.53), the result (6.14). 

7. Synthesis for Subordination 

The scaling relations a = Xr 13 for (rwi) and (rw 3 ) and a = jih a for (rw2) 
imply jih a = Xt 13 , the scaling relation (4.7) for well-scaled transition to the 
diffusion limit of (CTRW). Synthetizing the limiting processes y = y(t*) of 
(ru! 2 ) and t* = t*(t) of (ra 3 ) we get the limiting process x = x(t) of (CTRW) 
as x = y(t*(t)), which implies the integral formula of subordination for the 
solution of the Cauchy problem (2.1), namely 

poo 

u(x,t)= / f a ,e(x,Q q(t*,t) dU , (7.1) 
Ju=o 

with 

q(U,t) = q , Q (U,t) = t^tJ^L/ fcWt) = ' M P {U/1?) , (7.2) 

which is a probability density in t* evolving with t, see (6.13)-(6.14). Note 
that we denote the parent process by in distinction to the limit x(t) 

of (CTRW). Compare with Meerschaert et al. (2002) and with Gorenflo, 
Mainardi and Vivoli (2007). 
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The limit of (rwi) and (ra 2 ) yields our method of parametric subordination. 
They describe the processes t = t(t*) and y = y(t*), hence, by identifying 
y with the spatial variable x and eliminating i*, our method of parametric 
subordination x = y(t*), t = £(£*) for constructing trajectories of particles. 
In this method we call y = y(t*) the "parent process" , t = £(£*) the "leading 
process". Leaning on Feller (1971) we call t* = t*(t) the "directing process" . 

Let us remark that the methods described by Fogedby (1994), Kleinhans 
and Friedrich (2007), Zhang, Meerschaert and Baeumer (2008) working 
with Langevin equations are, in a certain sense, equivalent to parametric 
subordination. In the two coupled stochastic differential equations (6.3) 
and (6.4) of Gorenflo, Mainardi and Vivoli (2007) we have hinted to this 
Langevin approach without calling it so. Fogedby (1994) describes in detail 
this approach for the spatially multi-dimensional situation and discusses 
its relevance for treatment of problems from applications. Kleinhans and 
Friedrich (2007) produce and exhibit several particle trajectories and discuss 
information that can be drawn from these. Zhang, Meerschaert and Baeumer 
(2008) do likewise and consider furthermore situations of diffusion type and 
with variable coefficients in the spatial operator. As a concrete application 
they outline what is happening in Bromide plumes occurring in a site of North 
America, actual irregular processes that have inspired hydrological research. 
They compare results of measurements with those of their simulations and 
so demonstrate convincingly the power of their theory. 

It is our intention to compare our approach to simulation of particle 
trajectories with other approaches and generalizations in a forthcoming 
paper. The essence of our approach is that it is based on a systematic 
and consequent application of the CTRW integral equation to the various 
processes involved. 

8. Numerical Results 

In this Section, after describing the numerical schemes adopted, we shall 
show the sample paths for two case studies of symmetric (6 = 0) fractional 
diffusion processes: {a = 2, j3 = 0.80}, {a = 1.5, j3 = 0.90}. As explained 
in the previous sections, for each case we need to construct the sample paths 
for three distinct processes, the leading process t — t(t*) , the parent process 
x = y(t*) (both in the operational time) and, finally, the subordinated process 
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x = x(t), corresponding to the required fractional diffusion process. We shall 
depict the above sample paths in Figs. 1, 2, 3, respectively, devoting the left 
and the right plates to the different case studies. For this purpose, following 
Gorenflo, Mainardi and Vivoli (2007), we proceed as follows. 

First, let the operational time t* assume only integer values, say = n 
with n = 0, 1, ... , 10000. Then, produce 10000 independent identically 
distributed random deviates, say X±, X 2 , . . . , Xloooo, having a symmetric 
stable probability distribution of order a. Now, with the points 

n 

x = 0, x n = n>l, (8.1) 

k=l 

the couples (t* )n , x n ), plotted in the (t*,x) plane (operational time, physical 
space) can be considered as points of a sample path {?/(£*) : < t* < 10000} 
of a symmetric Levy motion with order a corresponding to the integer values 
of operational time = t* tn . In this identification of with n we use the 
fact that our stable laws for waiting times and jumps imply A = fi = 1 in the 
asymptotics (4.1) and (4.2) and r = h = 1 as initial scaling factors in (4.3). 

In order to complete the sample path we agree to connect every two successive 
points (£*,„,:£„) and (t^ n+ i,x n+ i) by a horizontal line from (£*,„, x„) to 
Xn), and a vertical line from (t^ n+ i,x n ) to (t*,„+i, x n+ i). Obviously, 
that is not the "true" Levy motion from point (£*,„, x n ) to point (t^ n+ i, x n+ i), 
but from the theory of CTRW we know this kind of sample path to converge 
to the corresponding Levy motion paths in the diffusion limit. However, as 
the successive values of t* i7l and x n are generated by successively adding the 
relevant standardized stable random deviates, the obtained sets of points 
in the three coordinate planes: (t*,t), (t*,x), (t,x) can, in view of infinite 
divisibility and self-similarity of the stable probability distributions, be 
considered as snapshots of the corresponding true random processes occurring 
in continuous operational time and physical time t, correspondingly. 
Clearly, fine details between successive points are missing. 

The well-scaled passage to the diffusion limit here consists simply in regularly 
subdividing the {t*} intervals of length 1 into smaller and smaller subintervals 
(all of equal length) and adjusting the random increments of t and x according 
to the requirement of self-similarity. Furthermore if we watch a sample path 
in a large interval of operational time the points (£*,„, x n ) and (t* 5 „ + i, x n+1 ) 
will in the graphs appear very near to each other in operational time 
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and aside from missing mutually cancelling jumps up and down (extremely 
near to each other) we have a good picture of the true processes. Plots in 
Fig. 1 (devoted to the leading process, the limit of (rwi)) thus represent 
sample paths in the (t*,t) plane of unilateral Levy motions of order (3. By 
interchanging the coordinate axes we can consider Fig 1 as representing 
sample paths of the directing process, the limit of (rw 3 ). 

Plots in Fig. 2 (devoted to the parent process, the limit of (rw 2 )) represent 
sample paths in the (t*,x) plane, produced in the way explained above, for 
Levy motions of order a and skewness 8 = (symmetric stable distributions). 
As indicated above, it is possible to produce 10000 independent identically 
distributed random deviates, say 7\, T 2 , . . . , Tioooo having a stable probability 
distribution with order j3 and skewness —f3 (extremal stable distributions). 
Then, consider the points 

n 

to = 0, t n = J2 T k, n>l, (8.2) 

k=l 

and plot the couples (i*, n ,i n ) i n the (£*,£) (operational time, physical time) 
plane. By connecting points with horizontal and vertical lines we get sample 
paths {t(t m ) : < t* < 10000} describing the evolution of the physical 
time t with the increasing of the operational time t*. Now, plotting points 
{t{t*, n ),y{t*,n)) i n the (t,x) plane, namely the physical time-space plane, and 
connecting them as before, one gets a good approximation of the sample 
paths of the subordinated fractional diffusion process of parameters a, (3 and 
= 0. 

In Fig. 3 (devoted to the subordinated process, the limit of (CTRW)) we 
show paths obtained in this way, from the points calculated in the previous 
paths. In Fig. 3 Left we plotted the points (i(i*, n ), y(t*, n )) obtained in Fig. 1 
Left and Fig. 2 Left, while Fig. 3 Right shows the points of Fig. 1 Right and 
Fig. 2 Right. 

By observing the figures the reader will note that horizontal segments 
(waiting times) in the (t, x) plane (Fig. 3) correspond to vertical segments 
(jumps) in the (t,t*) plane (Fig. 1). Actually, the graphs in the (t,x)- 
plane depict continuous time random walks with waiting times (shown as 
horizontal segments) and jumps Xk (shown as vertical segments). The left 
endpoints of the horizontal segments can be considered as snapshots of the 
true particle path (the true random process to be simulated), the segments 
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Figure 1: A sample path for (rwi), the leading process t = t(t*). 
LEFT: {a = 2, (3 = 0.80}, RIGHT: {a = 1.5 , = 0.90}. 




2000 4000 6000 8000 10000 2000 4000 6000 8000 10000 



t, t« 

Figure 2: A sample path for (ru^), the parent process x = y(t*). 
LEFT: {a = 2, (3 = 0.80}, RIGHT: {a = 1.5 , (3 = 0.90}. 




Figure 3: A sample path for (no), the subordinated process x = x(t). 
LEFT: {a = 2, (3 = 0.80}, RIGHT: {a = 1.5 , £ = 0.90}. 



being segments of our ignorance. In the interval t n < t < t n+1 the true 
process (namely the spatial variable x = X(t)) may jump up and down 
(infinitely) often, the sum (or integral) of all these ups (counted positive) 
and downs (counted negative) amounting to the vertical jump X n+ i. Finer 
details will become visible by choosing in the operational time i* the step 
length r << 1 (instead of length 1 as we have done) and correspondingly the 
waiting times and spatial jumps as r 1 ^ multiplied by a standard extreme /in- 
stable deviate, r 1//a multiplied by a standard (in our special case: symmetric) 
a-stable deviate, respectively, as required by the self-similarity properties of 
the stable probability distributions. In a forthcoming paper we will describe 
in detail what happens for finer and finer discretization of the operational 
time t*. 



Figure 4: Diagram for the connections between the four random walks (rwi), 
(rw 2 ),(rw 3 ) and (CTRW) = (rw). 

After sketching our method of well-scaled transition to the diffusion limit 
under power law regime, we have slpitted the CTRW into three random 
walks(rtyi), (™ 2 ), (rw 3 ), each of them involving two of the three directions: 
natural time, operational time and space. The random walk (rw 3 ) is the 
inverse of (ru^). To better visualize the matter, here Fig. 4 shows in form 
of a diagram the connections between the three random walks (rwi), (rtu 2 ), 
(rw 3 ) and (rw) = (CTRW). 



9. Conclusions 
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In section 8 the figures 1, 2, 3 have shown some numerical realizations of 
these random walks. We can consider Fig. 1 as a representation of (rwi), or 
by interchange of axes as one of (rw^), Fig. 2 as one of (r^), and finally 
Fig. 3 representation of (CTRW). 

We have carried out the analysis of transition by aid of the transforms of 
Fourier and Laplace, thereby consistently using corresponding variants of 
the Montroll- Weiss equation (not only in space and natural time, but also 
not conventionally in the above-mentioned combinations). In this way we 
have exhibited the true power of this equation in asymptotic analysis of 
random walks. In a quite natural way we have interpreted the diffusion limits 
as leading to the common subordination formula as well as the method of 
parametric subordination of constructing trajectories of particles. 
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